!>\brief
!! This is a simple test case on the unit square \f$[0\,,1]^d\f$.
!!
!! \n
!!
!! We consider the problem
!! \f{displaymath}{
!!  \begin{array}{rcll}
!!   \nabla\cdot(-\epsilon\nabla u + au) + \sigma u & = & f & {\rm
!!   in}\,\Omega \\\
!!   u&=&0 & {\rm on}\,\partial\Omega,
!!  \end{array}
!! \f}
!! with \f$\Omega=[0\,,1]^d\f$ and the coefficients are specified as
!! follows. Given the parameters \f$\epsilon\in\mathbb{R}^+\f$,
!! \f$m\in(\mathbb{N}^+)^d\f$ and the function
!! \f{displaymath}{
!!  \eta_s(t) = 1-e^{-\frac{1-t^s}{s\epsilon}}, \qquad s=1,2,\ldots
!! \f}
!! we let \f$a_i = x_i^{m_i-1}\f$, \f$\sigma = \sum_{i=1}^d
!! \sigma_i(x)\f$, with
!! \f{displaymath}{
!!  \sigma_i(x) = \left\{\begin{array}{ll}
!!   0, & {\rm if}\,m_i=1 \\\
!!   x_i^{m_i-2}, & {\rm if}\,m_i\geq2,
!!  \end{array}\right.
!! \f}
!! and \f$f = \sum_{i=1}^d f_i(x)\f$, with
!! \f{displaymath}{
!!  f_i(x) = \left\{\begin{array}{ll}
!!   \displaystyle
!!   (1+m_i-\eta_{m_i}(x_i))x_i^{m_i-1} \prod_{j\neq i} x_j\eta_{m_j}(x_j),
!!     & {\rm if}\,m_i=1 \\\
!!   \displaystyle
!!   (1+m_i)x_i^{m_i-1} \prod_{j\neq i} x_j\eta_{m_j}(x_j),
!!     & {\rm if}\,m_i\geq2.
!!  \end{array}\right.
!! \f}
!! The analytic solution of this test case is
!! \f{displaymath}{
!!  u(x) = \prod_{i=1}^d x_i\eta_{m_i}(x_i)
!! \f}
!! with gradient
!! \f{displaymath}{
!!  \partial_{x_i}u(x) = 
!!     \left( \eta_{m_i}(x_i)-\frac{x_i^{m_i}}{\epsilon}
!!      (1-\eta_{m_i}(x_i)) \right)\prod_{j\neq i} x_j\eta_{m_j}(x_j),
!! \f}
!! total flux
!! \f{displaymath}{
!!  q_i = \left(-\epsilon\eta_{m_i}(x_i)+x_i^{m_i}\right)
!!            \prod_{j\neq i} x_j\eta_{m_j}(x_j)
!! \f}
!! and divergence of the total flux
!! \f{displaymath}{
!!  \nabla\cdot q = \sum_{i=1}^d x_i^{m_i-1}(1+m_i-\eta_{m_1}(x_i))
!!     \prod_{j\neq i} x_j\eta_{m_j}(x_j).
!! \f}
!!
!! \note For \f$\epsilon\to0\f$ the analytic solution exhibits a sharp
!! boundary layer. Since the coefficients \f$\sigma\f$ and \f$f\f$
!! remain bounded for \f$\epsilon\to0\f$, the boundary layer is a
!! genuine hyperbolic effect.
!<----------------------------------------------------------------------
module mod_bubble_test

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_initialized, &
   new_file_unit

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_bubble_test_constructor, &
   mod_bubble_test_destructor,  &
   mod_bubble_test_initialized, &
   test_name, &
   test_description,&
   coeff_diff,&
   coeff_adv, &
   coeff_xiadv,&
   coeff_re,  &
   coeff_f,   &
   coeff_dir, &
   coeff_neu, &
   coeff_rob, &
   coeff_lam, &
   coeff_q,   &
   coeff_divq

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 character(len=*), parameter ::    &
   test_name = "bubble"

! Module variables

 ! public members
 character(len=100), protected ::    &
   test_description(3)
 logical, protected ::               &
   mod_bubble_test_initialized = .false.
 ! private members
 integer :: m(3) ! exponents
 real(wp) :: eps ! diffusion
 character(len=*), parameter :: &
   test_input_file_name = 'bubble_test.in'
 character(len=*), parameter :: &
   this_mod_name = 'mod_bubble_test'

 ! Input namelist
 namelist /input/ &
   m, eps

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_bubble_test_constructor(d)
  integer, intent(in) :: d

  integer :: fu, ierr
  character(len=9) :: des_format
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.)    .or. &
       (mod_fu_manager_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_bubble_test_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   ! Read input file
   call new_file_unit(fu,ierr)
   open(fu,file=trim(test_input_file_name), &
      status='old',action='read',form='formatted',iostat=ierr)
    if(ierr.ne.0) call error(this_sub_name,this_mod_name, &
      'Problems opening the input file')
    read(fu,input)
   close(fu,iostat=ierr)

   write(test_description(1),'(a,i1)') &
     'Bubble test case, dimension d = ',d
   write(des_format,'(a,i1,a)') '(a,',d,'i2,a)'
   write(test_description(2),des_format) &
     '  exponents m = [',m(1:d),' ];'
   write(test_description(3),'(a,e9.3,a)') &
     '  diffusion eps = ',eps,'.'

   mod_bubble_test_initialized = .true.
 end subroutine mod_bubble_test_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_bubble_test_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_bubble_test_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_bubble_test_initialized = .false.
 end subroutine mod_bubble_test_destructor

!-----------------------------------------------------------------------
 
 pure function coeff_diff(x) result(mu)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: mu(size(x,1),size(x,1),size(x,2))

  integer :: i
 
   mu = 0.0_wp
   do i=1,size(x,1)
     mu(i,i,:) = eps
   enddo
 end function coeff_diff
 
!-----------------------------------------------------------------------
 
 pure function coeff_adv(x) result(a)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: a(size(x,1),size(x,2))

  integer :: i

   do i=1,size(x,1)
     a(i,:) = x(i,:)**(m(i)-1)
   enddo
 end function coeff_adv
 
!-----------------------------------------------------------------------

 pure function coeff_xiadv(x) result(xi)
 ! Coefficient \xi such that
 !  a = \mu \nabla \xi
  real(wp), intent(in) :: x(:,:)
  real(wp) :: xi(size(x,2))

  integer :: i
  
   xi = 0.0_wp
   do i=1,size(x,1)
     xi = xi + (x(i,:)**m(i))/real(m(i),wp)
   enddo
   xi = xi/eps
 end function coeff_xiadv
 
!-----------------------------------------------------------------------
 
 pure function coeff_re(x) result(sigma)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: sigma(size(x,2))
 
  integer :: i

   sigma = 0.0_wp
   do i=1,size(x,1)
     if(m(i).ge.2) sigma = sigma + x(i,:)**(m(i)-2)
   enddo
 end function coeff_re
 
!-----------------------------------------------------------------------
 
 pure function coeff_f(x) result(f)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: f(size(x,2))

  integer :: i, imask(size(x,1),size(x,2))
  real(wp) :: p(size(x,1),size(x,2)), fi(size(x,2))

   do i=1,size(x,1)
     imask(i,:) = i
     p(i,:) = x(i,:)*eta(m(i),x(i,:))
   enddo
 
   f = 0.0_wp
   do i=1,size(x,1)
     if(m(i).eq.1) then
       fi = (real(1+m(i),wp) - eta(m(i),x(i,:)))*x(i,:)**(m(i)-1)
     else
       fi =  real(1+m(i),wp)                    *x(i,:)**(m(i)-1)
     endif
     
     f = f + fi*product(p,1,imask.ne.i)
   enddo
 end function coeff_f
 
!-----------------------------------------------------------------------

 pure function coeff_dir(x,breg) result(d)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: d(size(x,2))

   d = coeff_lam(x)
 
 end function coeff_dir
 
!-----------------------------------------------------------------------

 pure function coeff_neu(x,breg) result(h)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: h(size(x,2))

   h = 1.0_wp
 
 end function coeff_neu
 
!-----------------------------------------------------------------------

 pure function coeff_rob(x,breg) result(g)
  real(wp), intent(in) :: x(:,:)
  integer, intent(in) :: breg
  real(wp) :: g(size(x,2))

   g = 0.0_wp
 
 end function coeff_rob
 
!-----------------------------------------------------------------------

 pure function coeff_lam(x) result(l)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: l(size(x,2))

  integer :: i

   l = 1.0_wp
   do i=1,size(x,1)
     l = l * x(i,:)*eta(m(i),x(i,:))
   enddo
 end function coeff_lam
 
!-----------------------------------------------------------------------

 pure function coeff_q(x) result(q)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: q(size(x,1),size(x,2))

  integer :: i, imask(size(x,1),size(x,2))
  real(wp) :: p(size(x,1),size(x,2))

   do i=1,size(x,1)
     imask(i,:) = i
     p(i,:) = x(i,:)*eta(m(i),x(i,:))
   enddo
 
   do i=1,size(x,1)
     q(i,:) = (-eps*eta(m(i),x(i,:)) + x(i,:)**m(i)) &
                *product(p,1,imask.ne.i)
   enddo
 
 end function coeff_q
 
!-----------------------------------------------------------------------

 pure function coeff_divq(x) result(dq)
  real(wp), intent(in) :: x(:,:)
  real(wp) :: dq(size(x,2))

  integer :: i, imask(size(x,1),size(x,2))
  real(wp) :: p(size(x,1),size(x,2))

   do i=1,size(x,1)
     imask(i,:) = i
     p(i,:) = x(i,:)*eta(m(i),x(i,:))
   enddo
 
   dq = 0.0_wp
   do i=1,size(x,1)
     dq = dq + x(i,:)**(m(i)-1) * (real(m(i)+1,wp)-eta(m(i),x(i,:))) &
               * product(p,1,imask.ne.i)
   enddo
 end function coeff_divq
 
!-----------------------------------------------------------------------
 
 pure function eta(s,t)
  integer, intent(in) :: s
  real(wp), intent(in) :: t(:)
  real(wp) :: eta(size(t))
 
   eta = 1.0_wp - exp((t**s-1.0_wp)/(real(s,wp)*eps))
 end function eta
 
!-----------------------------------------------------------------------

end module mod_bubble_test

